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Abstract 

We present a dynamo mechanism arising from the presence of barotrop- 
ically unstable zonal jet currents in a rotating spherical shell. The 
shear instability of the zonal flow develops in the form of a global 
Rossby mode, whose azimuthal wavenumber depends on the width of 
the zonal jets. We obtain self-sustained magnetic fields at magnetic 
Reynolds numbers greater than 10 3 . We show that the propagation 
of the Rossby waves is crucial for dynamo action. The amplitude of 
the axisymmetric poloidal magnetic field depends on the wavenumber 
of the Rossby mode, and hence on the width of the zonal jets. We 
discuss the plausibility of this dynamo mechanism for generating the 
magnetic field of the giant planets. Our results suggest a possible link 
between the topology of the magnetic field and the profile of the zonal 
winds observed at the surface of the giant planets. For narrow Jupiter- 
like jets, the poloidal magnetic field is dominated by an axial dipole 
whereas for wide Neptune-like jets, the axisymmetric poloidal field is 
weak. 



1 Introduction 

The zonal (i.e. axisymmetric and azimuthally directed) jet streams visible at 
the surface of the giant planets are a persistent feature of the fluid dynam- 
ics of these planets (figure 1). The gas giants (Jupiter and Saturn) display 
a strong eastward equatorial jet, extending to latitudes ±20° with a peak 
velocity exceeding 100 m/s on Jupiter (Porco et al., 2003), and to latitudes 
±30° with a peak velocity exceeding 400 m/s on Saturn (Sanchez-Lavega 
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Figure 1: Zonal velocity measured at the surface in the planet's mean rotat- 
ing frame for each of the four giants by tracking cloud features in the outer 
weather layer. Profiles adapted from Porco et al. (2003), Sanchez-Lavega 
et al. (2000), Sromovsky et al. (2001) and Sromovsky and Fry (2005). 

et al., 2000). At higher latitudes, alternating prograde (eastward) and ret- 
rograde (westward) jets of smaller amplitude are observed extending all the 
way to the poles. These profiles are fairly symmetric with respect to the 
equator. On the ice giants (Uranus and Neptune) the picture is rather differ- 
ent. A very intense retrograde equatorial current is present with maximum 
velocity of 100 m/s on Uranus (Sromovsky and Fry, 2005) and 400 m/s on 
Neptune (Sromovsky et al., 2001). At higher latitudes, a single prograde 
jet of large amplitude is present in each hemisphere. Several decades of ob- 
servations show that these zonal flows remain approximately steady (Porco 
et al., 2003). 

The origin of these zonal flows and the associated question of the depth 
to which they extend into the planets' interiors have been areas of ac- 
tive research in rotating fluid dynamics for several decades (e.g. Jones and 
Kuzanyan, 2009, and references therein; see also the review by Vasavada and 
Showman, 2005). In particular, several models have been proposed to ex- 
plain the zonal wind pattern of Jupiter, and can be categorized into two main 
classes: weather layer models and deep convective layer models. The former 
assume that the zonal flows are produced in a shallow stably stratified region 
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near cloud level. These models are able to reproduce the high latitude struc- 
tures with alternating eastward and westward jets and a strong equatorial 
current (e.g. Williams, 1978; Cho and Polvani, 1996). These models tend 
to produce a retrograde equatorial jet (Yano et al., 2003), so they provide a 
plausible explanation for the retrograde equatorial flow of the ice giants but 
not for the prograde flow observed on gas giants. A parametrized forcing 
such as a strong equatorially-localized baroclinicity is required to force a 
shallow system to produce a prograde equatorial jet (Williams, 2003). The 
second class of models is deep convection models which simulate most or 
all of the whole 10 4 km-thick molecular hydrogen layer (Busse, 1976; Chris- 
tensen, 2001, 2002; Manneville and Olson, 1996). The presence of deep 
convection is inferred from the observation that the atmospheres of the ma- 
jor planets emit more energy by long-wave radiation than they absorb from 
the Sun. Consequently their atmospheres must receive additional heat sup- 
plied by the interior of the planet. Recent numerical models using either 
a Boussinesq approximation (Heimpel et al., 2005) or an anelastic approx- 
imation (Jones and Kuzanyan, 2009) and low Ekman numbers (i.e. strong 
rotational effect compared with viscous dissipation) display alternating zonal 
jets at high latitudes. A strong eastward equatorial jet is a robust feature 
of these models where the Coriolis force dominates buoyancy, in good agree- 
ment with the gas giant observations. Interestingly, deep convection models 
suggest that the zonal velocity generated by non-linear interactions of con- 
vective motions (i.e. the motions directly forced by buoyancy) is roughly 
geostrophic, that is, invariant along the direction of the rotation axis. This 
feature is also present in strongly compressible models provided that the 
Ekman number is small enough, despite the increase of density with depth 
yielding ageostrophic convective motions (Jones and Kuzanyan, 2009; Kaspi 
et al., 2009). When the convection is more vigorous such that the buoyancy 
force overcomes the Coriolis force, 3D turbulence homogenizes angular mo- 
mentum; a retrograde jet forms in the equatorial region and a single strong 
prograde jet forms in the polar region, in good agreement with the ice giant 
observations (Aurnou et al., 2007). 

Another feature of the giant planets is their strong magnetic fields (fig- 
ure 2). The observed magnetic fields for gas and ice giants differ drastically 
(see for instance the recent review by Russell and Dougherty, 2010). Jupiter 
and Saturn have a main axial dipole component (corresponding to Z = 1, 
m = in figure 2), a feature shared with the Earth for instance (Yu et al., 
2010; Burton et al., 2009). Neptune and Uranus, on the other hand, have 
strong non-axial multipolar components (corresponding to I = 2, 3 in fig- 
ure 2) compared with the axial dipole component (Connerney et al., 1991; 
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Figure 2: Spectra of the magnetic field squared amplitude at the planetary 
radius for degrees I and order m up to 3 obtained from inversion models 
of the magnetic measurements. The squared amplitude for a given degree 
I is A\ = Y^m=Q^ + 1) [fe" 1 ) 2 + (hf) 2 ] using a Schmidt normalisation for 
the spherical harmonics. The squared amplitude for a given mode m is 

A m = El=m( Z + l ) [iaD 2 + ( h F) 2 ]- 9™ and h T are the Gauss coefficients 
in gauss. After Yu et al. (2010) (Model Galileo 15), Burton et al. (2009) 
(Cassini measurements), Connerney et al. (1991) (model Os) and Herbert 
(2009) (AH5 model from magnetic observations and auroral data). 
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Herbert, 2009). The magnetic field is generated in the deep, electrically 
conducting regions of the planets' interiors: a metallic hydrogen layer for 
Jupiter and Saturn (Nellis et al., 1999; Guillot, 2005, and references therein) 
and an electrolyte layer composed of water, methane and ammonia (Hub- 
bard et al., 1991; Nellis et al., 1997) or superionic water (Redmer et al., 
2011) for Uranus and Neptune. 

Numerical models of convective dynamos in rapidly rotating spherical 
shells typically produce axial dipolar dominated magnetic fields for moder- 
ate Rayleigh numbers and moderate Ekman numbers (e.g. Olson et al., 1999; 
Aubert and Wicht, 2004; Christensen and Wicht, 2007). To explain the un- 
usual large scale non-dipolar magnetic fields of Uranus and Neptune, models 
using peculiar parameter regimes or different convective region geometries 
have been proposed. The latter models show that a numerical dynamo op- 
erating in a thin shell surrounding a stably-stratified fluid interior produces 
magnetic field morphologies similar to those of Uranus and Neptune (Hub- 
bard et al., 1995; Holme and Bloxham, 1996; Stanley and Bloxham, 2006). 
Gomez-Perez and Heimpel (2007) obtain weakly dipolar and strongly tilted 
dynamo magnetic fields when high magnetic diffusivities are used (or equiv- 
alently small electrical conductivity). Their results show that these peculiar 
fields are stable in the presence of strong zonal circulation and when the flow 
has a dominant effect over the magnetic fields. This feature is also empha- 
sized by Aubert and Wicht (2004) who find stable equatorial dipole solutions 
with a weak magnetic field strength and low Elsasser number (measure of the 
relative importance of the Lorentz and Coriolis forces) for moderately low 
Ekman numbers. They argue that the magnetic field geometry of the equa- 
torial dipole solution is incompatible with the columnar convective motions 
and thus this morphology is stable only when Lorentz forces are weak. 

Although scaling laws derived from numerical simulations of dynamos 
driven by basal heating convection predict dipolar magnetic field in plan- 
etary parameter regimes (Olson and Christensen, 2006), recent numerical 
simulations using more realistic parameter values (lower Ekman numbers) 
have not produced large scale magnetic fields so far, and require larger mag- 
netic Reynolds numbers (measure of magnetic induction versus magnetic 
diffusion) (Kageyama et al., 2008). Moreover, convection in the interior of 
Jupiter is often thought to be driven by secular cooling (Stevenson, 2003). 
Numerical dynamos driven by secular cooling typically produce weak dipole 
or multipolar magnetic field for larger forcing (Kutzner and Christensen, 
2000; Olson and Christensen, 2006) depending on boundary conditions (Hori 
et al., 2010). Therefore the question of the generation of large scale magnetic 
field by turbulent convective motions in the planetary parameter regime re- 
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mains open. 

The dichotomies observed in the magnetic fields and in the zonal wind 
profiles of the giant planets are rather striking. Up to now no study has tried 
to relate them directly, probably because the former is a feature of the deep 
interior whereas the latter is a characteristic of the surface. However, if some 
mechanism is able to transport angular momentum from the surface down 
to the deep, fully conducting region then the zonal motions may influence 
the generation of the magnetic field. In the non-magnetic deep convection 
models (Heimpel et al., 2005; Jones and Kuzanyan, 2009), zonal motions 
extend geostrophically throughout the electrically insulating molecular hy- 
drogen layer down to the bottom of the model. On the other hand, due 
to the possible rapid increase of electrical conductivity with depth in the 
outer region, Liu et al. (2008) argued that the ohmic dissipation produced 
by geostrophic zonal motions shearing dipolar magnetic field lines would ex- 
ceed the luminosity measured at the surface of Jupiter if the vertical extent 
of this geostrophic zonal motions exceeds 4% of the planet radius. However, 
the argument of Liu et al. (2008) is purely kinematic, that is the action of 
the magnetic forces on the flow and the feedback on the magnetic field are 
ignored. In a self-consistent magnetohydrodynamic model, the zonal flow 
would adjust toward a non- geostrophic state due to the action of magnetic 
forces if the electrical conductivity of the fluid is significant (Glatzmaier 
(2008), see also the non-linear numerical simulations of convectively-driven 
dynamos of Aubert (2005)). In this case, angular momentum may be trans- 
ported along the magnetic field lines leading to a dynamical state close to 
the Ferraro state. This state minimizes the ohmic dissipation produced by 
the shearing of the poloidal magnetic field by the zonal flow as the poloidal 
magnetic field lines are aligned with angular velocity contours. Both scenar- 
ios, either geostrophic zonal balance or Ferraro state, imply the existence of 
multiple zonal jets of significant amplitude at the top of the fully conducting 
region beneath. The plausibility of each scenario depends on the radial pro- 
file of electrical conductivity, which is currently not well constrained within 
the giant planets (Nellis et al., 1999). 

The idea of the work presented in this paper is that these zonal jets may 
exert, by viscous or electromagnetic coupling, an external forcing at the top 
of the deeper conducting envelope. From previous studies (Schaeffer and 
Cardin, 2006; Guervilly and Cardin, 2010) we know that the viscous cou- 
pling between a differentially rotating boundary and a low-viscosity electri- 
cally conducting fluid can generate a self-sustained magnetic field in different 
geometries. Zonal motions can be subject to barotropic shear instabilities 
which have a lengthscale independent of the viscosity, unlike convective in- 
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stabilities. These instabilities are able to generate large scale magnetic fields, 
and so they are an interesting source of dynamo action under planetary in- 
terior conditions. In order to test the plausibility of a dynamo driven by this 
source in isolation, we use an incompressible 3D numerical dynamo model 
with a zonal velocity profile imposed at the top of a spherical shell contain- 
ing a conducting fluid. We use a dynamical approach, that is non-linear 
interactions between the flow and the magnetic field are taken into account; 
therefore the fluid flow is free to adopt a three-dimensional structure as long 
as it satisfies the imposed viscous boundary conditions. 

The dynamics of the deep conducting region is usually assumed to be 
slower than the dynamics of the outer molecular hydrogen region due to 
magnetic braking, even if uncertainties remain in the electrical conductivity. 
The model presented in this paper assumes an idealized one-way coupling 
between the outer and deep regions. A more realistic model would need to 
account for the back reaction of the deep layer onto the outer layer; a study 
of the consistent dynamical interaction of the two layers is beyond the scope 
of this paper. For studies of more realistic coupling, see promising recent 
numerical models of self-consistent convectively-driven dynamos in spheri- 
cal shells including radially variable electrical conductivity of Heimpel and 
Gomez Perez (2011) and Stanley and Glatzmaier (2010). In these models, 
slow convective motions in the interior dynamo region coexist with strong 
zonal flow near the outer surface. Differential rotation in the interior is only 
partially inhibited by the strong magnetic field. 

In order to assess the role of the zonal wind profile on the topology of 
the sustained magnetic field, we use both Jupiter-like and Neptune-like zonal 
wind profiles. In the giant planets, as in rocky planets, it is usually assumed 
that the dynamo mechanism is driven by convective motions. The giant 
planets display a strong surface heat flux (with the exception of Uranus) 
meaning that heat transfer is efficient in the interior of the planet and thus 
mostly due to convection (Guillot and Gautier, 2007, and references therein). 
Here we want to assess the efficiency of zonal velocity forcing alone, so we 
do not model convective motions. 

The first goal of this work is to quantify what amplitude of the zonal wind 
inside the conducting layer is needed to trigger the dynamo instability, so we 
do not model the exact or realistic coupling between the molecular hydrogen 
upper layer and the deep, electrically conducting region. Our second goal is 
to test to what extent the pattern of the zonal flow imposed at the top of 
the conducting layer influences the topology of the self-sustained magnetic 
field. 

We first describe the model and the numerical method used (section 2). 
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Then we present numerical results from simulations in the non-magnetic case 
(section 3) followed by results from dynamo simulations (section 4). The 
application of our results to planetary conditions is discussed in section 5. 

2 Model 

We model the deep conducting layer of the giant planets as a thick spher- 
ical shell. At the top of the conducting layer we impose an axisymmetric 
azimuthal velocity to represent the zonal flow generated in the overlying en- 
velope. The shell rotates around the z-axis at the imposed rotation rate fi. 
The aspect ratio is 7 = rj/r Q where rj is the inner sphere radius, correspond- 
ing to a rocky core, and r Q the outer sphere radius, corresponding to the 
top of the fully conducting region. The fluid is assumed incompressible with 
constant density p and constant temperature, that is, no convective motions 
are computed. The assumption of incompressibility is made for simplicity, 
although the pressure scale height at the depths of the conducting layer is 
roughly 8000km (Guillot et al., 2004), that is, about 1/5 of the thickness of 
the layer. The effects of compressibility may well play a role in the dynamics 
of the conducting regions (see for instance Evonuk and Glatzmaier, 2004). 

For simplicity we model the angular momentum coupling with the ex- 
ternal zonal flow as a rigid boundary condition for the velocity at the outer 
boundary, rather than as a shear stress condition. The flow is driven through 
a boundary forcing rather than a volume forcing to avoid directly impos- 
ing bidimensionality to the velocity field. As we are interested in the bulk 
magnetohydro dynamical process, the exact nature of the coupling (electro- 
magnetic or viscous, shear stress or rigid) with the upper molecular hydrogen 
layer is not crucial for our study. We discuss the implication of the choice 
of the rigid boundary condition in section 3. The radial profile of electrical 
conductivity is not well constrained in the gas giants. In particular the ex- 
istence of a first order or continuous transition between the molecular and 
metallic hydrogen phase is still an open question, although high-pressure 
experiments are in favor of a continuous transition (Nellis et al., 1999). We 
choose to model the outer boundary as electrically insulating to simplify the 
coupling between the layers. The conductivity is assumed constant through- 
out the whole modeled conducting layer. As we do not model the molecular 
hydrogen layer, we assume zonal geostrophic balance within this envelope 
for simplicity. The amplitude of the zonal motions at the outer boundary of 
our model is therefore the same as the surface winds. This idealized repre- 
sentation of the dynamics of the molecular hydrogen layer would be altered 
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if the magnetic forces upset the zonal geostrophic balance. Depending on 
the magnitude and radial profile of the electrical conductivity, the amplitude 
of the zonal motions might be reduced, and the zonal flow contours would 
tend to align with the magnetic field lines, although we do not expect the 
characteristics of the zonal jets (narrow or wide, relative amplitude of the 
peaks) to be altered very much. 

We use two different synthetic azimuthal velocity profiles for the bound- 
ary forcing imposed at the top: a multiple jet profile for the gas giants with 
a profile based on Jupiter's surface zonal winds (hereafter profile J) and a 
3-band profile based on Neptune's surface zonal winds (profile N). 

For Jupiter, we use the profile given in Wicht et al. (2002) 



where s = r sin 9, r s is the surface radius of the planet and Uq = U(ro, 9 = vr/2) 
no controls the numbers of jets. The profile at the radius r s best matches 
the observed profile at the surface for uq = 4 (figure 3). The profile U(r Q , 9) 
is used to drive the flow at the top of our simulated metallic hydrogen layer 
(figure 3). The ratio j s = r s /r Q determines the U profile at r Q . We choose 
7s = r s /r = 1/0.8 = 1.25 following Guillot et al. (1994). 

For the Neptune-like profile, we use the zonal velocity profile measured 
at the surface of Neptune, approximated by a polynomial of order 10 in 
latitude. We project this surface velocity profile geostrophically down to r Q 
using 7s = 1/0.85 = 1.18 (Hubbard et al, 1991) (figure 3). 

The existence of a rocky core at the centre of the giant planets is un- 
certain and depends on the poorly constrained composition of the planet. 
Estimates for the core mass are — 14m® for Jupiter (total mass 318771®), 
6— 17m® for Saturn (total mass 95m®) and 0— 4m® for Uranus and Neptune 
(total mass 15m,® and 17m® respectively) where m® denotes the mass of 
the Earth (Guillot, 2005). If present, the rocky cores are therefore believed 
to be small. Following the interior model of Jupiter proposed by Guillot 
et al. (1994) we use an aspect ratio rj/r D = 0.2 for all the simulations per- 
formed. The inner core is assumed to be electrically conducting, with the 
same conductivity as the fluid in the conducting layer. We did not carry out 
simulations with an insulating core as the effect of the conductivity of the 
inner core on the dynamo mechanism is believed to be small (Wicht, 2002). 
The velocity boundary condition is no-slip at the inner boundary. 

The velocity u is scaled by Uq, the absolute value of the azimuthal ve- 
locity imposed at the equator of the outer sphere. The lengthscale is the 
radius of the outer sphere r Q . The magnetic field B is scaled by \J pnor ClUo 
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Figure 3: Zonal velocity profile imposed at the surface of model J (left) and 
model N (right) (solid lines). Both profiles are obtained by assuming that 
the zonal velocities are geostrophic for r s > r > r Q and using the profile 
represented by a dashed line at the surface of the planet (r = r s ): model 
J, profile (1) with no = 4, 7 S = r s /r Q = 1.25 and Uq = 100; model N: 
polynomial fit of order 10 in latitude of the zonal wind profile measured at 
the surface of Neptune (figure 1) with 7 S = 1/0.85 = 1.18. For comparison 
the zonal wind profile measured at the surface of Jupiter is plotted in gray. 
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where p is the fluid density and fj>o is the vacuum magnetic permeability. 
We numerically solve the momentum equation for an incompressible fluid, 



Re^ + Re (u • V) u + -e z x u = -Vp + V 2 u + — (V x B) x B, (2) 

ot hi h 

the continuity equation, 

V • u = 0, (3) 

and the magnetic induction equation, 

i9B 1 

V • B = 0, (5) 

where p is the dimensionless pressure, which includes the centrifugal poten- 
tial. 

The Reynolds number Re = r Uo/v parametrizes the mechanical forcing 
exerted on the system by controlling the amplitude of the zonal velocity. 
The magnetic Prandtl number Pm = u/rj measures the ratio of viscous 
to magnetic diffusivities. The magnetic Reynolds number Rm is defined 
as Rm = RePm. The Ekman number E = u/(Qr^) measures the im- 
portance of the viscous term over the Coriolis force. The Rossby number 
Ro = ReE = Uo/(Qr Q ) is the ratio of inertial force to Coriolis force. Note 
that in our definition the Rossby number refers to the amplitude of the pre- 
scribed zonal jets at the surface, and not to the local flow velocity. 
The results presented in this paper were obtained with the PARODY code, 
a fully three-dimensional and non-linear code. The code was derived from 
Dormy (1997) by J. Aubert, P. Cardin, E. Dormy in the dynamo benchmark 
(Christensen et al., 2001), and parallelised and optimised by J. Aubert and 
E. Dormy. The velocity and magnetic fields are decomposed into poloidal 
and toroidal scalars and expanded in spherical harmonic functions in the 
angular coordinates with I representing the latitudinal degree and m the 
azimuthal order. A finite difference scheme is used on an irregular radial 
grid (finer near the boundaries to resolve the boundary layers). A Crank- 
Nicolson scheme is implemented for the time integration of the diffusion 
terms and an Adams-Bashforth procedure is used for the other terms. 
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3 Dynamics without the magnetic field 



For a rapidly rotating system in which the Coriolis force exactly balances 
the pressure force, the Proudman- Taylor constraint states that the flow is 
z-invariant and follows geostrophic contours. For an incompressible fluid 
in a bounded container, these geostrophic contours correspond to surfaces 
of equal height. In a sphere the only geostrophic motions are azimuthal 
and axisymmetric. In the giant planets' conducting envelopes, the Ekman 
number is about 10 -16 and the Rossby number is much smaller than 1 
(Guillot et al., 2004). In the absence of a magnetic field, we expect the 
Proudman- Taylor constraint to hold for large scale motions. As we want 
to reach the dynamical regime in which the flow is strongly geostrophic, 
the use of small Ekman and Rossby numbers is required. We carried out 
simulations for 10~ 5 > E > 10~ 6 for model J and 10~ 5 >£>5x 10~ 6 for 
model N. The Rossby numbers are always smaller than 0.1. For the profile 
J, in cases of low Ekman numbers (E < 2 x 10~ 6 ), we imposed longitudinal 
symmetry by calculating only the harmonics of a chosen order m s . The 
required resolution for E = 10 -6 is 500 points on the radial grid and I = 580 
spherical harmonics degrees. 

3.1 Axisymmetric flow 

When the imposed boundary forcing is small enough, i.e. when the Rossby 
number Ro is less than a critical value Ro c , the flow is axisymmetric and 
predominantly azimuthal (figure 4). The zonal jets imposed at the outer 
boundary extend into the volume along lines parallel to the axis of rotation. 

The use of no-slip boundary conditions yields a differential rotation be- 
tween the boundary and the bulk of the fluid. This differential rotation 
is accommodated across viscous Ekman boundary layers, which scale as 
(E/ cos^) 1 / 2 , where 9 is the colatitude. By Ekman pumping, viscous forces 
within the Ekman layers drive axial motions of order E 1 ! 2 within the bulk 
of the fluid (figure 4). These meridional circulations advect angular momen- 
tum from the boundary layer into the bulk of the fluid and cause the jets to 
propagate faster than by pure viscous diffusion. At low latitudes, the Ekman 
layer is thicker so the Ekman pumping is stronger, yielding to a more effi- 
cient driving of the zonal motions in the bulk by the outer boundary layer. 
For model J (figure 5(a)), the zonal velocity in the bulk relative to that 
imposed at the outer boundary is noticeably weaker for the inner jets than 
for the outer jets. When E decreases this effect is less marked, and in the 
E — > limit we expect the basic zonal velocity to be perfectly geostrophic 
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Figure 4: Angular velocity u^j '(r s'mO) (left) and streamlines of the merid- 
ional circulation (isocontours of tjj = rsmO^f- with u p the velocity poloidal 
scalar) (right) of the axisymmetric flow in the northern meridional plane. 
For the meridional circulation, anti-clockwise (clockwise) flows are shown in 
solid (dotted) lines. The parameter for the simulations are E = 5 x 1CP 6 
and Ro = 0.015 for model J (a) and E = 10~ 5 and Ro = 0.02 for model N 
(b). 
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Figure 5: Zonal velocity in the equatorial plane for subcritical numerical 
simulations (solid lines) compared to the zonal velocity at radius r = 0.98r o 
(symbols) and imposed velocity at the top (dashed line) both projected in 
the equatorial plane for (a) model J (E = 5 x 10 -6 (bold solid line and open 
squares) and E = 10~ 6 (thin solid line and black circles)) and (b) model N 
(E = 10~ 5 (bold solid line and open squares)). 

in the whole volume. The comparison between the zonal velocity just below 
the Ekman layer and in the equatorial plane (figure 5) shows that the zonal 
velocity is geostrophic in the bulk of the fluid (outside of the boundary lay- 
ers). For model N (figure 5(b)), the zonal jets are wider, so the zonal flow 
already displays a strong geostrophic structure at E = 10~ 5 . Note that the 
azimuthal velocity has to match the no-slip boundary condition at the inner 
core, and so an internal Stewartson layer forms on the axial cylinder tangent 
to the inner core (Stewartson, 1966). 

3.2 Non-axisymmetric motions 
3.2.1 Model J 

Rossby wave at the onset When the boundary forcing (measured by 
Ro) becomes greater than a critical value Ro c , the axisymmetric basic flow 
becomes unstable to a non-axisymmetric shear instability. The saturated 
instability takes the form of an azimuthal necklace of cyclonic and anticy- 
clonic vortices aligned with the axis of rotation, is nearly ^-independent and 
drifts eastward (figure 6). Close to the threshold, the radial extension of the 
pattern is large and occupies almost half of the gap. The pattern drifts with 
the same speed over its whole radial extension, even though the advection 
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Figure 6: Non-zonal axial vorticity in the equatorial plane (right) and in a 
meridional slice (left) for model J at E = 4 x 10 -6 and Ro = 1.01i?o c (blue: 
negative and red: positive). The black curve represents the zonal velocity 
in the equatorial plane. 

by the zonal flow velocity varies with s, implying that it is a single wave. 

Wicht et al. (2002) studied the linear stability of the imposed zonal 
flow (1) in a spherical shell modeling the insulating molecular hydrogen layer 
of Jupiter (aspect ratio 0.8). For E = 10 -4 they found nearly bidimensional 
instabilities that they described as drifting columns aligned with the rotation 
axis and similar to convective solutions. Although they do not identify 
these instabilities as waves, their characteristics are very similar to the ones 
obtained with our non-linear model. 

The nearly z-invariant structure and the prograde drift are two character- 
istics of Rossby waves propagating in a spherical container. The dispersion 
relation for the Rossby wave given by a local linear analysis is (e.g. Finlay, 



where (3 = h~ l (dh/ds) = —s/(r 2 — s 2 ) is related to the slope of the upper 
boundary of the spherical container of height h. k s and m/s are the ra- 
dial and azimuthal wavenumbers respectively. The theoretical Rossby wave 
frequency ui rw can be calculated at a given radius assuming k s « m/s and 
using the wavenumber m obtained from the numerical simulation. For dif- 
ferent E, the frequency to of the propagating wave observed in our numerical 
simulations always falls in the range cu rw (si) < u < uj rw (s2) where s\ = 0.56 



2008) 




(6) 
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(s2 = 0.87) is the smallest (resp. largest) radius where a significant vortic- 
ity associated with the presence of the wave can be seen in the numerical 
calculations. This strongly indicates that the shear instability occurs as a 
Rossby wave. 

The velocity of the zonal flow U enters the dispersion relation of the 
Rossby wave through a Doppler shift 

u(s) = u rw (s) + U(s) — . (7) 
s 

As reported earlier, oo(s) is constant in our numerical calculations so uj rw (s) 
must adapt in the s-direction for the wave to be coherent. In a prograde 
jet U > 0, uj rw must decrease, which requires a local increase in k s in 
equation (6) and so a local decrease in the radial lengthscale, which can 
be observed in figure 6. For small enough Ekman number (in practice 
E < 5 x 10~ 6 ), the critical wavenumber m c of the Rossby mode is inde- 
pendent of E. The radial lengthscale is determined by the width of the 
jet and the vortices are roughly circular in the equatorial plane (figure 6) 
suggesting that m c is controlled by the width of the jets. 

In a local approximation that neglects the curvature terms, a criterion of 
instability of barotropic shear flows has been derived by Ingersoll and Pollard 
(1982) for an anelastic model in a full rotating sphere and by Kuo (1949) for 
thin stably stratified "weather" layers. Using an inviscid Boussinesq model 
and for barotropic instability of a zonal flow U in a sphere, this necessary 
condition implies a change of sign of a quantity A at some radius: 

A = 2(3- Ro^, (8) 
ds 

where £ is the vorticity of the zonal flow, 

+ (9, 
as s 

Note that the curvature terms have been taken into account here. In a 
sphere, (3 is negative. Consequently, the zonal velocity profile is more prone 
to instability where the gradient of zonal vorticity is maximum and negative. 
Then for a profile U of sinusoidal form, the first shear instability occurs at the 
maximum of the prograde jets, and thus, perhaps surprisingly, at a null value 
of the zonal velocity shear dU /ds. Note that our numerical simulations show 
instabilities with a large radial extent and with maximum amplitude located 
in a retrograde zonal jet (see figure 6), even though the local instability 
criterion predicts an onset in a prograde jet. This observation emphasizes 
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that the local criterion does not predict the location of global saturated 
modes. 

The theoretical critical Rossby number obtained from applying the crite- 
rion (8) to the profile (1) imposed at the top of model is Ro^! 1 = 0.0011. The 
threshold of the first instability of the axisymmetric flow, denoted Ro™ hn , 
obtained with the numerical simulations are shown in figure 7. Despite the 
decrease of Ro™ hn with the Ekman number, Ro^ hn is still about four times 
larger than Ro^ for E = 10~ 6 because the amplitude of the zonal flow within 
the bulk is reduced by viscous boundary layers in the numerical simulations. 
Due to computational limitations, it is not possible for us to carry out simu- 
lations at smaller E with a fully non-linear code and prove the existence of an 
asymptotic regime for the inviscid instability threshold. For this purpose we 
used a dedicated linear code described in A. The linear code calculates lin- 
ear perturbation solutions to the momentum equation using the geostrophic 
profile U as the basic flow in the bulk of the fluid. The computational time 
is greatly reduced by the linear approach but is restricted to an analysis 
of the instability threshold. The growing solutions obtained with the linear 
code exhibit very similar features to the Rossby waves in the non-linear sim- 
ulations (frequency, bidimensional structure, radial extent, location of the 
maximum amplitude in a retrograde jet). In figure 7 the threshold Ro l * n 
obtained with the linear code approaches asymptotically the value given by 
the local theory. For the same Ekman number, Ro l ^ n is smaller than Ro™ lm 
since the geostrophic zonal flow U is used in the linear code, that is the jets 
in the bulk have greater amplitude than in the non-linear code. From our 
linear computations we conclude that the theoretical criterion (8) is relevant 
to explain the onset of instability obtained numerically. More details about 
the onset of the hydrodynamic instability can be found in Guervilly (2010). 

The characteristic time of the Rossby wave is r rw = 1/oj. At the instabil- 
ity threshold, the numerical simulations g 180" 1 for E < 5 x 10" 6 . 
The timescale of the zonal jets is T z j = r /Uo = il" 1 / Ro. For Ro = 0.01, we 
have r z j > T rw : the Rossby wave propagation is faster than the advection 
of the fluid by the zonal flow. The turnover time of a fluid particle trapped 
in a Rossby wave is Tt = l/V s where / is the typical radial displacement 
of the particle and V s the typical cylindrical radial velocity of the particle. 
At Ro = 1.01i?o c , V s is typically lO" 2 ^. In a rough approximation we use 
/ = 5, where 5 is the width of the jets, 5 ~ 0.1r o for the profile J. Then we 
obtain r to « 0.1r o /(10~ 2 C/ ) « WRo' 1 ^ 1 ps 10 3 ft -1 : the turnover time of 
the particle is much longer than the timescale of the wave. Consequently 
the particle oscillates rapidly as the wave propagates and is slowly advected 
by the zonal flow. In practice the radial displacement I is typically smaller 
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Figure 7: Critical Rossby number obtained from fully non-linear numerical 
simulations for model J (Ro^ hn , circles) compared to the theoretical Rossby 
number obtained with the local instability criterion (8) using the geostrophic 
profile (1) U(s,6 = vr/2) (-Ro*' 1 , black line). The critical Rossby number 
obtained from the linear numerical calculation is also shown (Ro^ 11 , crosses). 

than 5 and so the turnover time is slightly overestimated here. 

Supercritical regime When the Rossby number is increased in the su- 
percritical regime, other prograde jets will eventually become unstable. A 
second Rossby wave appears in the weakly supercritical regime, at Ro = 
1.06i?o c for E = 5 x 10 -6 , with a maximum velocity located in the retro- 
grade zonal jet at larger radius than the first wave maxima (i.e. the wave 
appearing for Ro = Ro c ) (figure 8(a)). To fill the larger circumference at 
larger radius the instability has a slightly larger wave number, m = 22 in- 
stead of 21, while the radial width of the jet is comparable. The second 
wave propagates faster, in agreement with the Rossby wave dispersion re- 
lation (6). Barotropic instabilities tend to broaden and weaken narrow jets 
by redistributing potential vorticity (see for instance Pedlosky, 1979). The 
smoothing of the jets saturates the amplitude of the Rossby waves. For this 
slightly supercritical regime the zonal flow profile is only weakly modified. 
Upon further increasing the forcing (Ro = 2.94Ro c ), several Rossby waves of 
different wavenumbers superpose and interact (figure 8(b)). The structure 
of the waves and the jets is still mainly bidimensional except in the viscous 
boundary layers. The typical cylindrical radial velocity is V s ~ O.lf/o and 
the Rossby number is 0.05 so the turnover time is about 20f2~ 1 assuming 
that the radial displacement I = 5, about the same order of magnitude as 
the timescale of the zonal jets. 

In figure 9(a) the time-averaged zonal flow in the equatorial plane is 
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(a) Ro = 1.06.Roc 




-0.1 -0.05 0.05 0.1 -0.2 -0.1 0.1 0.2 

(b) Ro = 2.94Ro c 



Figure 8: Snapshots of the radial (left) and azimuthal (right) velocity com- 
ponents in the equatorial plane for E = 5 x 10~ 6 and Ro > Ro c for model 
J. The velocities are scaled by Uq. For the colorscale has been truncated 
(u<p(r ,6 = 7r/2, 0) = 1). The black curve represents the zonal velocity in 
the equatorial plane. 
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(a) (b) 



Figure 9: (a) Time-averaged zonal velocity in the equatorial plane for 
E = 5 x 10~ 6 and different forcings, (b) Amplitude of the non-axisymmetric 
radial velocity V s (squares), non-axisymmetric azimuthal velocity V<f> (cir- 
cles), non-axisymmetric velocity (K? + V^) 1 / 2 (diamonds) and zonal veloc- 
ity at the radius s = 0.75 (triangles). All velocities were measured in the 
equatorial plane in the units of Uq. The amplitude of the non-axisymmetric 
velocity corresponds to the maximum in a snapshot, whereas the zonal ve- 
locity amplitude has been averaged in time. 

plotted for different Ro up to Ro = 5.88Ro c . As the forcing is increased, 
the Rossby waves gradually reduce the jet strength and broaden the jet 
width. For Ro = 2.94i?o c , the retrograde jet at s = 0.81 has been mostly 
destroyed leading to the widening of the zonal jet width. We note that 
the zonal flow becomes mostly westward for the strongest forcings. The 
amplitude of the zonal flow located at s > 0.9 is hardly affected because 
the threshold to destabilise the outermost jets is high due to the large slope 
(related to [3 in equation (8)). For Ro < 2.35i?o c , the amplitude of the non- 
axisymmetric velocity, relative to Uq, increases with the forcing (figure 9(b)). 
After reaching a maximum, at Ro = 2.35i?o c , the amplitude of the non- 
axisymmetric flow decreases relative to Uq. The "efficiency" of the forcing 
to drive the non-zonal velocity is reduced as the Rossby waves smooth the 
gradient of vorticity and so affect their excitation mechanism. 

The back reaction on the forcing velocity in the upper molecular hy- 
drogen layer is not taken into account in our model although it might sig- 
nificantly affect the zonal profile in the upper layer in the case of strong 
forcing. 
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Figure 10: Non-zonal axial vorticity in the equatorial plane (right) and in a 
meridional slice (left) for model N at E = 5 x 10 -6 and Ro = 1.01i?o c (blue: 
negative and red: positive). The black curve represents the zonal velocity 
in the equatorial plane. 

3.2.2 Model N 

The shear instability takes the form of an m = 2 oscillation in the azimuthal 
direction (figure 10). It is a single wave propagating eastward with the 
same frequency over the shell, and is nearly z-invariant. The maxima of 
the non-zonal vorticity are located on each side of the prograde jet. The 
characteristics of this wave are similar to the Rossby wave obtained with 
model J. The frequency of this wave is in agreement with the frequency of 
a theoretical Rossby wave of wavenumber m = 2 propagating at a radius 
s = 0.53 (assuming that k s ~ m/s in the dispersion relation (6)). For 
E = 10~ 5 and E = 5 x 10 -6 , the critical Rossby numbers obtained with 
the non-linear numerical simulations are respectively Ro™ hn = 0.0335 and 
R nhn _ q Q325. Using the instability criterion (8) with the profile imposed 
at the surface we obtain a critical Rossby number of 0.026 in good agreement 
with the non-linear numerical results when the Ekman number decreases. 

4 Magnetic field generation 

The non-axisymmetric motions are of prime importance for the dynamo 
mechanism because a purely toroidal flow cannot generate a self-sustained 
magnetic field. We note that some axisymmetric poloidal flow is present 
when Ro < Ro c as a weak meridional circulation is created by the Ekman 
pumping. However these axisymmetric motions are weak at small Ekman 
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numbers so we do not expect to find dynamos when the zonal flow is stable, 
that is when Ro < Ro c , in the asymptotic inviscid regime. Indeed we did 
not find dynamos when Ro < Ro c (up to Pm = 10). The non-axisymmetry 
associated with the hydrodynamic shear instability is a crucial element for 
the dynamo process: the stable zonal flow cannot sustain a magnetic field 
by itself. This is in agreement with the results obtained by Guervilly and 
Cardin (2010) with dynamos generated by spherical Couette flows (differen- 
tial rotation between two concentric spheres). 

4.1 Characteristics of the magnetic field for model J 

We have performed dynamo simulations for Ro = 1.17 — 1.76Ro c and E = 
5 x 10~ 6 . We find that the dynamo threshold occurs at a rather high value 
of the magnetic Prandtl number, Pm c ~ 5. The critical magnetic Reynolds 
number (defined via the maximum forcing velocity) required for dynamo 
action is Rm c ~ 20, 000 (see section 5 for an estimate of the critical magnetic 
Reynolds number defined via the local velocity). For a given forcing, we have 
performed calculations just above the critical magnetic Prandtl number, 
Pm c , and up to 2Pm c . 

The main features of the self-sustained magnetic field can be observed 
in figures 11 and 12. The magnetic field displays a dipolar symmetry, i.e. 
antisymmetry with respect to the equatorial plane, 

(B r ,B e ,B^(r,7r-e,^) = (—B r , Bg, —B^r, 9, 4>). (10) 

The magnetic field is predominantly toroidal and axisymmetric (correspond- 
ing to the mode m = in figure 12(a)). The toroidal magnetic field does 
not emerge from the conducting region as the outer region is electrically 
insulating. The strongest poloidal component is the axial dipole within the 
conducting region and outside of the outer sphere (corresponding to the 
harmonic (l,m) = (1,0) in figure 12). Within the bulk of the flow, the ax- 
isymmetric poloidal magnetic field lines are mostly significantly bent where 
the Rossby wave causes a strong magnetic induction (figure 11(a)). A mag- 
netic field at the scale of the Rossby wave is produced in this region as can 
be observed on the spectra of magnetic energy (figure 12(a)) with signifi- 
cant peaks at m = 22 in the poloidal and toroidal magnetic energies and 
at I = 23 in the poloidal magnetic energy (I — m is odd to preserve the 
dipolar symmetry). Close to the outer boundary, the axisymmetric poloidal 
magnetic field lines converge and diverge locally (figure 11(a)). This is due 
to the induction of axisymmetric magnetic field by the secondary meridional 
circulation produced by Ekman pumping. This effect is very localized and 
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(a) Axisymmetric magnetic field 




(b) Radial magnetic field at r s 

Figure 11: Magnetic field for model J. (a) Snapshot of the axisymmetric 
magnetic field in a meridional plane: magnetic poloidal field lines (left) and 
azimuthal magnetic field (right) (blue: negative and red: positive), (b) 
Map of the radial magnetic field at the surface of the planet r s = 1.25r Q in 
unit of 10~ 3 y / p/Io^ r o (solid line: positive and dotted line: negative). The 
poloidal magnetic field at r = r s is calculated assuming the region between 
r Q and r s is electrically insulating. The parameters of this simulation are 
E = 5x 10" 6 , Ro = 1.17Ro c and Pm = 5« Pm c . 
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(a) Kinetic and magnetic energies in the fluid 




(b) Ai and A m at r s 



Figure 12: Magnetic energy spectra for model J. (a): Kinetic (KE) and 
magnetic (ME) energy per unit volume for each spherical harmonics degree 
I (left) and mode m (right) in the fluid conducting region given in unit 
of pUq. (b): Squared amplitudes of the magnetic field, A[ (left) and A m 
(right) as defined in figure 2, at r s = 1.25r D given in unit of ppoU^. Only 
the degrees of significant amplitude have been plotted, that is, I even for 
the kinetic poloidal and magnetic toroidal energies and I odd for the kinetic 
toroidal and magnetic poloidal energies. These data are taken at a particular 
instant and have not been time-averaged. Same parameters than figure 11. 
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generates a magnetic field of small latitudinal scale that decreases rapidly 
with radius. 

The spectrum and map of the radial magnetic field at the surface of our 
modeled planet (at radius r s = 1.25r G ) (Figs. 11(b) and 12(b)) show that 
the magnetic field is strongly dominated by the axial dipole. The magnetic 
field generated at the scale of the Rossby wave (m = 22) is still visible 
in the spectrum of the magnetic field but its amplitude is weak at this 
radius: about four orders of magnitude smaller than the amplitude of the 
axisymmetric mode (note that the spectrum in figure 12(b) represents the 
squared amplitude of the field). 

In all the simulations performed, no inversion of polarity of the axial 
dipole has been observed. The tilt of the dipole is rather weak, at most 2° 
from the rotation axis. We found a secular variation of the dipole axis of 
about 1° every 1000 rotation periods or alternatively 0.001 global magnetic 
diffusion time. 

Just above the dynamo threshold (Pm c < Pm ^ 2Pm c ), the magnetic 
field is weak: the magnetic energy contained within the fluid conducting 
region is only about 5% of the kinetic energy. The magnetic field does 
not strongly act back on the flow, except to produce its own saturation. 
A comparison between the zonal flow in the non-magnetic case and in the 
presence of the dynamo magnetic field does not reveal significant differences. 
The magnetic field lines of the poloidal field are almost aligned with the 
rotation axis and the flow structure (see figure 11(a)) so the flow disruption 
due to Lorentz forces is weak. 

4.2 Characteristics of the magnetic field for model N 

We performed simulations at Ro = 1.05 — 1.5i?o c and E = 10 -5 . We find 
the dynamo threshold at Pm c ~ 1, that is, the critical magnetic Reynolds 
number is Rm c ~ 4000. 

The main features of the self-sustained magnetic field can be observed in 
figures 13 and 14. The self-sustained magnetic field displays an equatorial 
symmetry, i.e. 

(B r , B e , B^){r, tt - 0, <£) = (B r , —B e , B^)(r, 9, <j>). (11) 

Within the fluid conducting region, the axisymmetric toroidal field is the 
strongest component whereas the poloidal field is dominated by the m = 2 
mode, not the axisymmetric m = mode. The m = 2 mode corresponds to 
the magnetic field generated at the scale of the Rossby wave (figure 14(a)). 
The axisymmetric poloidal field is multipolar, mainly composed by the 
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(a) Axisymmetric magnetic field 




(b) Radial magnetic field at r s 

Figure 13: Magnetic field for model N (same as figure 11). For the ax- 
isymmetric azimuthal field, blue corresponds to negative values and red to 
zero values. The radial magnetic field is plotted at the surface of the planet 
r s = 1.18r Q in unit of 10 -5 y / p/Io£^o- The parameters of this simulation are 
E = 1(T 5 , Ro = 1.20Ro c and Pm = 2 » 2Pm c . 
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(a) Kinetic and magnetic energies in the fluid 




I m 

(b) Ai and A m at r s 



Figure 14: Magnetic energy spectra for model N (same as figure 12). Only 
the degrees of significant amplitude have been plotted, that is, I even for the 
kinetic poloidal and magnetic poloidal energies (plus I = 1) and I odd for 
the kinetic toroidal and magnetic toroidal energies. Same parameters than 
figure 13. 
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(l,m) = (2,0) (axial quadrupole) and (l,m) = (4,0) modes. At the sur- 
face of the planet (figure 14(b)), the magnetic field appears to be mainly 
axisymmetric (with the I = 2 and I = 4 harmonics degrees dominant). The 
amplitude of the m = 2 structure is weak, about two orders of magnitude 
smaller than the m = mode but still visible at high latitudes on the map 
of the radial field at the surface (figure 13(b)). 

Due to the equatorial symmetry of the field, the magnetic field lines in 
the equatorial plane are roughly perpendicular to the cylindrical structure 
of the flow whereas they are nearly aligned at higher latitudes (figure 13(a)). 
As a result magnetic braking acting on the flow is stronger in the equato- 
rial region than at high latitude regions. In the simulation performed here 
(Pm c < Pm ^ 2Pm c ), the magnetic energy is weak compared to the ki- 
netic energy (about 5%) so the feedback of the magnetic field on the flow 
remains weak. For a stronger magnetic field (at larger magnetic Reynolds 
numbers), we expect that the flow disruption would become important. As 
a result the equatorially symmetric solution may become unstable and the 
magnetic field may switch to an axial dipolar symmetry. This is the result 
obtained by Aubert and Wicht (2004) in convectively-driven dynamos: they 
found equatorial dipolar magnetic fields for Rayleigh numbers close to the 
convection onset; these solutions become unstable as the convective forcing 
is increased and an axial dipolar configuration is preferred. 

In summary, the flows driven by the profiles J and N produce very dif- 
ferent poloidal magnetic fields: mainly a strongly axisymmetric dipole for 
the profile J and a weak multipolar axisymmetric field dominated by the 
magnetic field induced at the scale of the Rossby waves for the profile N. 
In both cases the magnetic field within the conducting region is mainly 
an axisymmetric toroidal field. The different magnetic field morphology is 
quite surprising given that the flows are quite similar: strong zonal flows 
and propagating Rossby waves. In the next section we review the dynamo 
mechanism that has been proposed to operate for similar flows and suggest 
the key difference between profiles J and N that determines the topology of 
their self-sustained magnetic fields. 

4.3 Dynamo mechanism 

Using a quasi-geostrophic flow and a kinematic approach (no Lorentz force in 
the momentum equation), Schaeffer and Cardin (2006, hereafter SC06) ob- 
tain numerical dynamos generated by an unstable axisymmetric shear layer 
(Stewartson layer): for a strong enough forcing, the Stewartson layer is un- 
stable to non- axisymmetric shear instabilities, which appear in the form of 



28 



Rossby waves (of wavenumber about 10 for the Ekman numbers and Rossby 
numbers they investigated). The self-sustained magnetic field has a strong 
axisymmetric toroidal component and a mostly axisymmetric poloidal com- 
ponent. SC06 show that the time dependence of the flow is a key ingredient 
for the dynamo effect: time-stepping the magnetic induction equation using 
a steady flow taken either from a snapshot or a time-average leads to the 
decay of the magnetic field. They characterize the dynamo process as an 
aw mechanism. In mean field theory, the a effect parameterizes the gen- 
eration of an axisymmetric poloidal magnetic field from the correlation of 
small scale magnetic field and velocity. The a effect usually requires that 
the flow possess some helicity, the correlation between fluid velocity and 
vorticity, H = u u (e.g. Moffatt, 1978). Flows displaying a columnar struc- 
ture aligned with the axis of rotation, such as Rossby waves or convection 
columns (Olson et al., 1999), typically possess strong mean helicity. As 
these columns are essentially bidimensional vortical structures, the helicity 
is mainly produced by the term u z uj z . In nearly z-invariant flow, the axial 
(z) velocity is mostly due to two terms: the slope effect and the Ekman 
pumping. The slope effect comes from the combination of mass conserva- 
tion and impenetrable boundaries: a (cylindrical) radial velocity u s creates 
an axial velocity u z ~ zf3u s with j3 = h~ 1 (dh/ds). In the limit of rapid 
rotation in a spherical container, this contribution is much larger (of order 
1) than the Ekman pumping (u z ~ E 1 I 2 uj z ). However, the axial velocity 
produced by the slope effect is phase shifted by ir/2 with uj z , and so does 
not allow the production of mean helicity. On the contrary, axial velocity 
produced by Ekman pumping is in phase with the axial vorticity and a dy- 
namo mechanism based on the Ekman pumping associated to an azimuthal 
necklace of axial vortices is plausible (Busse, 1975). In a numerical experi- 
ment at small Ekman numbers (E = 0(1O~ 8 )), SC06 artificially remove the 
Ekman pumping and observe dynamo action with nearly the same threshold 
showing that the Ekman pumping is unimportant in their dynamo mecha- 
nism. The crucial importance of the time dependence of the flow and the 
negligible contribution of the Ekman pumping lead SC06 to consider the in- 
volvement of the Rossby waves in the dynamo process. They conjecture that 
the propagation of the Rossby waves yields a proper phase shift between the 
non-axisymmetric magnetic field and velocity field in order to produce the 
axisymmetric poloidal magnetic field. 

Avalos-Zuniga et al. (2009) have calculated the a tensor, describing the 
generation of a large scale magnetic field by correlation of small scale velocity 
and magnetic field, with a flow geometry corresponding to Rossby waves. In 
the absence of Ekman pumping, they show that the diagonal components of 
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the a tensor, which are the relevant coefficients for the a effect, are non-zero 
if and only if the flow pattern is drifting relative to the mean flow. 

Tilgner (2008) explains that the time dependence of a velocity field can 
lead to dynamo action even when any particular snapshot of the velocity field 
cannot because the linear operator associated with the induction equation 
is non-normal. In particular, he shows that the simple time dependence of 
a propagating wave is enough for dynamo action. Several numerical studies 
report the importance of the time dependence of the velocity field, mainly 
of oscillating nature (Reuter et al., 2009; Gubbins, 2008). 

The idea that the propagation of Rossby waves may maintain a dynamo 
action is very appealing as their presence is ubiquitous in rotating fluid dy- 
namics. A system in which no wave propagation occurs, and which is unable 
to produce u z by another mechanism, such as buoyancy, will rely on Ekman 
pumping to create axial velocity with the proper phase shift. However, in 
the limit of small Ekman number, the Ekman pumping vanishes and the 
dynamo threshold should become infinitely high. The dynamo mechanism 
relying on the propagation of Rossby waves is robust in the limit of small 
Ekman number as the presence of these waves does not rely on the action 
of viscosity. 

Due to the close resemblance of the flow (zonal motions and propa- 
gating Rossby wave) in our 3D numerical model and the kinematic quasi- 
geostrophic model of SC06, we now try to establish if the dynamo mechanism 
evoked in SC06 is at work in our 3D model. To formalize their idea, let us 
first consider a simple theoretical model. The velocity field is composed by 
a zonal flow, U(s)ej,, and the small scale velocity of a Rossby wave u m with 

u m (s, </>, z, t) = (u™(s, z)e s + u$(8, z)e + uf(s, z)e z )e i ^-^ (12) 

where u™, u™ and u™ are complex and uj is the frequency of the wave. 
The magnetic field is composed of an axisymmetric magnetic field B, and 
a magnetic field perturbation induced at the scale of the Rossby wave b m 
with 

b m (s, 0, z, t) = (bf(s, z)e s + b%(s, z)e^ + b™(s, z)e z )e i ^-^ +xt (13) 

where b™ , 6™ and b™ are complex and A is the growth rate of the magnetic 
field. The equations for the evolution of the poloidal components of B in 
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cylindrical coordinates B s and B z are 



dt 



d ( B \ 

— (ufbf - u™b™) + r, (V 2 B~ S - -|J , (14) 



dB v 1 d 



__ s - + V V 2 B Z , (15) 

where the overbar denotes an azimuthal average. It is immediately apparent 
that if u™ (u™) is out of phase by ir/2 with b™ (resp. 6™), then B s and B z 
will be decaying in time. If we suppose that B^ 3> B S ,B Z the equations for 
b™ and 6™ are 

(A-ic-)&™ = — u^ + ? ? V 2 6™. (17) 
s s 

where c = (u/(m/s) — U) is the phase speed of the wave relative to the 
mean flow U. In the case of marginal stability (A = 0), if we neglect the 
magnetic diffusivity r\ then we obtain that b™ (b™) is in phase with u™ 
(u™ resp.). Moreover if the axial velocity is mainly due to the slope effect 
then u™ = zfiu™ and so according to the equations (16)-(17) 6™ ~ zflb™. 
This implies that the first term of the right hand side of equations (14)-(15) 
is almost zero and thus B s and B z are decaying. Consequently magnetic 
diffusivity at the scale of b™ and 6™ must play a role in the generation of 
the axisymmetric poloidal magnetic field by introducing a short phase lag 
between the velocity and magnetic modes. This phase lag depends on the 
spatial structures of b™ and and hence the terms u^b™ and uV^b™ do 
not cancel out. Note that the importance of magnetic diffusivity is well 
established in the a effect (Roberts, 2007). On the other hand, if the wave 
is not propagating, c = 0, then 



im 



[ ? db r L l /)™ s 



-™<*B* = >)V 2 K". (19) 

In this case the magnetic field perturbations b™ and b™ are out of phase 
with u™ and u™ (as u™ and u™ are correlated by the slope effect) and so 
the averaged products u r z a b T ^ and u™b™ are zero. We can conclude that in 
order for this simple model to work as a mean-field dynamo (i) the wave 
must propagate and (ii) the magnetic diffusivity must act on the magnetic 
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(a) Model J (b) Model N 



Figure 15: Non-axisymmetric magnetic field (coloured) and non- 
axisymmetric velocity (black lines: positive, and gray lines: negative) in 
a plane a few degree of latitude above the equatorial plane (northern hemi- 
sphere). Same parameters than figures 11 and 13. 

field generated at the scale of the waves. As the Rossby wave propagates, 
the location of the induction of the magnetic field perturbation is forced to 
drift with the same rate, but with a phase-shift. The phase-shift between 
the magnetic field perturbation and the Rossby wave depends on both the 
phase speed c and the magnetic diffusivity rj. The argument above implies 
that this phase-shift is essential for the dynamo mechanism. 

Using any particular snapshot of the velocity field for time stepping the 
magnetic induction in our numerical simulations with models J or N leads 
to the decay of the magnetic field. The failure of dynamo in the kinematic 
numerical experiment with both models is readily explained by our simple 
theoretical model. 

In figure 15, we plot the non-axisymmetric components of the velocity, 
u™ and u™ and magnetic field, b™ and b™ obtained in the numerical simula- 
tions for model J and model N in a plane located just above the equatorial 
plane (b™ and b™ are zero in the equatorial plane by dipolar symmetry in 
model J). The correlation of u™ with u™ confirms that u™ is mainly pro- 
duced by the slope effect for both models. For model J (figure 15(a)), we 
observe that and 6™ are in phase so has a significant amplitude. 

However, b™ is out of phase with u™, which means that u™b™ 3> u™b™. 
This may be an effect of the magnetic diffusivity as b™ and b™ have dif- 
ferent spatial structures, or due to radial derivatives of B s and B z that 
we neglect in equation (17). Consequently u™b™ mainly contributes to the 
generation of strong B s and B z . 

For model N (figure 15(b)) strong positive (negative) crescent-shaped 
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patches of b™ and b™ are visible in the cyclonic (resp. anticy clonic) vor- 
tices, out of phase by ir/2 with u™ and u™. Consequently these crescent- 
shaped structures of b™ and 6™ do not contribute to the terms u™b™ and 
ymfom ^he presence of these maxima of b™ and b™ are not explained by the 
theoretical model (equations (16) and (17)) likely because of the neglect of 
the axial and radial derivatives of B s and B z , which are important in this 
region (see figure 13(a)). Round-shaped lobes of b™ and b" 1 of weaker ampli- 
tude (located in the middle of the gap) are observed in phase with u™ and 
u™. Consequently these round-shaped structures of b™ and b™ contribute 
to the terms u™b™ and u™b™. Unlike model J (where u™6™ 3> u™b™), 
u™b™ ~ u™6™ so only a weak axisymmetric multipolar magnetic field is 
maintained in this case. At the surface of the planet this axisymmetric field 
is the dominant component but in comparison with the strongly axisymmet- 
ric dipolar field produced in model J, the field is of small amplitude: the 
amplitude of the axisymmetric radial field is about 10~ 3 y/p[ioUo for model 
J at Rm = 1. 17 Rm c (Ro = 1.17 Ro c and Pm « Pm c ) (figure 12(b)) while 
it is only 10" 5 ^/ppqUq for model N at Rm c = 2ARm c (Ro = 1. 20 Ro c and 
Pm = 2Pm c ) (figure 14(b)). 

The main difference between the Rossby waves in models J and N is their 
size. The phase speed of the Rossby wave, c « tt(3/(m/s) 2 , is about 100 
times larger for a m = 2 wave than a m = 22 wave, for a fixed radius s and 
rotation rate f2. On the other hand, the magnetic diffusion acts more rapidly 
on small scale structures. The typical propagation timescale for a Rossby 
wave of size d is r rw = 1/(0, /3d) assuming that the radial and azimuthal 
lengthscales of the wave are similar. The magnetic diffusion timescale at 
the scale of the vortex d is t v = d 2 /rj. The ratio of the two timescales is 

= m. (20) 
Trw V 

The dependence to the third power of the size, d oc 1/m, shows that the 
magnetic diffusion timescale relative to the propagation timescale is about 
three orders of magnitude smaller for an m = 22 mode than anra = 2 mode 
for the same parameter values. For the simulation presented for model N, the 
ratio Tr)/T rw is about 10 5 . For model J the ratio t^/t^ is about 500 so the 
propagation of the Rossby wave is still much more rapid than the magnetic 
diffusion. For both models, we found that the values of the small scale 
magnetic field in phase with the velocity is of the same order of magnitude. 
The velocity field of the vortices is also about the same order of magnitude 
for the two models. The difference between the two models is that, in model 
N, the magnetic diffusion acts too slowly on the m = 2 magnetic structures 
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compared to the wave propagation to produce a significant enough phase 
lag between 6™ (b™) and u™ (u™). Consequently, the term u™b™ — uV}b™ is 
weak and leads to little generation of axisymmetric poloidal magnetic field. 

The last stage of the dynamo mechanism is the generation of the ax- 
isymmetric toroidal field. It can either be produced from the correlation of 
small scale velocity and magnetic field (as an a effect) or an uj effect, that is 
the shearing of the axisymmetric poloidal magnetic field by the mean zonal 
flow U. SC06 find that the uj effect from the Stewartson layer is domi- 
nant in their numerical model. The zonal shear produced in the Stewartson 
layer is stronger than the shear we obtained with the profiles J and N, so 
it is not clear that the uj effect is important in our model prima facie. In 
a 2 dynamos, both toroidal and poloidal components are typically of sim- 
ilar magnitudes (Olson et al., 1999). Here, the strong toroidal magnetic 
field suggests that the uj effect is more important. To confirm this, we plot 
in figure 16 the term responsible for the uj effect in the azimuthal com- 
ponent of the magnetic induction equation (Gubbins and Roberts, 1987), 
rB r d r {r~ l U) + r~ 1 sm#f?0<90(sin# _1 i7). For model J, as we expect, this 
term is most significant in the region where the poloidal magnetic field lines 
are bent and misaligned with the zonal flow structure (see figure 11(a)). The 
correlation of sign and location of the maxima of the uj effect in the bulk of 
the fluid with the axisymmetric azimuthal field indicates that it is mainly 
generated by the uj effect. Note that some uj effect is also present close to 
the outer boundary, where the poloidal magnetic field lines converge and 
diverge locally due to induction by the Ekman pumping. However no par- 
ticularly strong axisymmetric azimuthal magnetic field is produced in this 
region (figure 11(a)) so this small scale field diffuses probably very rapidly. 
For model N the outer part of the jet (s > 0.5) is retrograde and creates a 
negative uj effect whose sign and location correlate with the axisymmetric 
azimuthal field, implying that the main dynamo process in the outer region 
is indeed the u effect. However, and the uj effect are anti-correlated in 
the inner region (s < 0.5) so another dynamo process such as a correlation 
of small scale velocity and magnetic field must be at work there. 

We have not yet addressed the question of the selection of the axial 
dipolar symmetry or the axial quadrupolar symmetry. In kinematic dy- 
namo calculations, Gubbins et al. (2000) show that minor changes in the 
flow can select very different eigenvectors. For a self-consistent system the 
selection rules are thus very subtle. As discussed in section 4.2, Aubert 
and Wicht (2004) found that axial quadrupolar symmetry is incompatible 
with the vertical structures of cyclones and anticyclones in convectively- 
driven dynamos, and so these solutions are unstable for strong convective 
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(a) Model J (b) Model N 



Figure 16: uo effect in the meridional plane in the bulk (outside the Ekman 
layers) (blue: negative and red: positive). Same parameters than figures 11 
and 13. 

flows. In our simulations of model N, this conclusion suggests that the ax- 
ial quadrupolar symmetry would be unstable for larger magnetic Reynolds 
numbers, and an axial dipolar field would be preferred. The selection of a 
given symmetry does not modify our argument that the wavenumber of the 
Rossby mode determines the amplitude of the axisymmetric magnetic field 
since no particular latitudinal symmetry is assumed. 

In this study, it appears that the dynamo mechanism relies on a subtle 
balance between the Rossby wave propagation and the magnetic diffusion 
and therefore is closely related to the size of the Rossby waves. The dynamo 
field produced with this mechanism requires high magnetic Reynolds num- 
bers {Rm c « 20,000 for model J and Rm c ~ 4000 for model N). However, 
in the limit of small Ekman number, this dynamo mechanism is expected 
to keep a finite value of the critical magnetic Reynolds number (Schaeffer 
and Cardin, 2006), whereas for dynamos that rely on Ekman pumping the 
critical magnetic Reynolds number becomes infinitely high. 

5 Summary and discussion 

We have numerically studied the dynamics of zonal flows driven by differen- 
tial rotation imposed at the top of a conducting layer and how they sustain 
a magnetic field. 
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5.1 Hydrodynamical instability 

In our hydrodynamical simulations, we found that the destabilisation of the 
zonal flow takes the form of a global (large radial extension) Rossby mode, 
even though the instability threshold is governed by a local criterion. The 
wavenumber depends on the width of the jets, and is independent of the 
viscosity and rotation rate provided that the former is sufficiently small. 
In the supercritical regime, several Rossby waves appear and saturate the 
amplitude of the zonal flow in the bulk of the fluid. They produce a widening 
of the jets and a strong damping of their amplitude, even for relatively small 
supercritical forcing (Ro = 2.94i?o c ). 

5.2 Constraints on the dynamo mechanism 

In the limit of small Ekman number, we find that the Rossby wave ap- 
pears for Ro c ~ 0.001 for a Jupiter-like zonal wind profile (model J) and 
Ro c ~ 0.02 for a Neptune-like profile (model N). In our numerical calcula- 
tions, non-axisymmetric motions are necessary for dynamo action to occur. 
As the viscosity is large in the numerical simulations compared to the plan- 
etary values, the Reynolds number is much smaller in the simulations. To 
reach a sufficiently high magnetic Reynolds number, the magnetic Prandtl 
number is of order 1, much larger than the expected planetary values. The 
critical magnetic Reynolds number Rm c is about 20, 000 for model J and 
4000 for model N. To make this dynamo mechanism work, two constraints 
must be satisfied: (i) Ro > Ro c and (ii) Rm > Rm c . Equivalently this 
gives constraints on the amplitude of the zonal motions at the top of the 
conducting region, Uq > Ro c Qr D , and on the electrical conductivity within 
the conducting region, a > Rm c /(Uor /j,o). 

The extrapolation of the constraint (i) to the giant planets is straight- 
forward as the hydrodynamical instability threshold is independent of the 
Ekman number, which is of order 10 -15 — 10~ 16 for Jupiter (Guillot et al., 
2004) and Neptune (Stevenson, 1983). For Jupiter (r G 56,000 km and 
fi = 1.8 x 10~ 4 s~ 1 ), the equatorial velocity at the top of the conducting 
region, Uq, must be larger than 10 m/s to have Ro > Ro c = 0.001. For 
Neptune (r 21,000 km and SI = 1.08 x 10 _4 s _1 ), Uq must be larger than 
45 m/s to have Ro > Ro c = 0.02. For both cases, this constraint is quite 
strong as it only allows for a factor 10 decrease of the amplitude of the zonal 
wind between the surface of the planet and the top of the deep conducting 
region, independently of the location of the top of this region. 

The extrapolation of the constraint (ii) to the giant planets requires 
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knowing how the critical Reynolds number scales with the Ekman number. 
When varying the Ekman number from 10 -6 down to 10 -8 , Schaeffer and 
Cardin (2006) found that Rm c remains constant (of the order of 10 4 in their 
simulations, close to the values found in our study). Based on their results, 
we assume that Rm c is of the same order of magnitude when the Ekman 
number is close to the planetary values. For Jupiter, we obtain that the 
electrical conductivity should be larger than 30S/m to have Rm > Rm c = 
20, 000 (using Uq = 10 m/s). For Neptune, the electrical conductivity should 
be larger than lOS/m to have Rm > Rm c = 4000 (using Uq = 45 m/s). 
This constraint on the conductivity is less restrictive than the constraint 
on the amplitude of the zonal motions and should be satisfied in the deep 
conducting layer of Jupiter (Nellis et al., 1999) and Neptune (Nellis et al., 
1997). 

We conclude that the differential rotation imposed by the zonal winds 
at the top of the conducting regions is a plausible candidate to drive the 
dynamo mechanism in the giant planets although a strong constraint on 
the amplitude of the zonal jet applies. Given the assumptions used in our 
model, such as incompressibility, constant conductivity, unrealistically large 
viscosity and viscous coupling between electrically insulating and conducting 
regions, this conclusion remains tentative. However, the robust nature of 
Rossby waves in the asymptotic limit of small Ekman numbers makes this 
dynamo mechanism appealing for planetary physical conditions. 

5.3 Generation of the axisymmetric field and width of the 
jets 

With a simple theoretical model, we show that the production of the ax- 
isymmetric field depends on the propagation of the Rossby waves and on 
the magnetic diffusion acting at the scale of the vortices. This model is in 
agreement with our numerical results: the magnetic diffusion rate of the 
m = 2 magnetic structures induced by the Rossby waves in model N is 
nearly negligible compared to the propagation rate of the wave: as a re- 
sult a weak axisymmetric poloidal magnetic field is generated; the magnetic 
diffusion acting on the m = 22 magnetic structures is not negligible com- 
pared to the propagation rate of the small size (m = 22) Rossby wave of 
model J: a dominant axisymmetric poloidal magnetic field is therefore gen- 
erated. Consequently, in this model, the width of the zonal jets has an 
important influence on the generation of the axisymmetric magnetic field 
by controlling the size of the Rossby waves. Our results suggest that the 
difference in the magnetic fields and the surface zonal winds may be related 
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if a (hydrodynamic or magnetohydrodynamic) mechanism can transport an- 
gular momentum between the surface and the deep, electrically conducting 
region. 

The critical magnetic Reynolds number of this dynamo mechanism is 
large. However, in order to compare with other dynamos, a more significant 
number may be the critical local magnetic Reynolds number associated with 
magnetic induction by the Rossby wave velocity Rm l c = V s djr\ where V s 
is the typical non-axisymmetric radial velocity and d is the lengthscale of 
the Rossby mode. For the dynamo obtained in model J (Rm c = 20,000), 
V s ~ O.lC/o an d m = 22 so we find Rm l c ~ 570. For the dynamo obtained in 
model N {Rm c = 4000), V s « 0MU and m = 2 so Rm l c » 130. Thus Rm l c 
is roughly 2—10 times larger than the magnetic Reynolds number needed for 
dynamo action with a convective forcing (Christensen and Aubert, 2006). 

5.4 Magnetic field at the planets' surfaces 

In our numerical model, we obtain a peak at small azimuthal scale in the 
magnetic field spectrum correlated with the width of the hydro dynamically 
unstable zonal jets. This is a testable prediction as the magnetic mea- 
surements of the forthcoming Juno mission (arrival at Jupiter in 2016) are 
expected to be of extraordinary quality due to the absence of a crustal mag- 
netic field on Jupiter. 

For model J, we obtain a secular variation of the dipole tilt of about 1° in 
1000 rotation periods or equivalently 0.001 global magnetic diffusion time. 
The dipole is strongly axisymmetric with a tilt that does not exceed 2°. On 
Jupiter, the dipole axis tilt measured with the Pioneer and Voyager data 
compared with the Galileo measurements is larger (about 10°) and displays 
a secular variation of about 0.5° in 20 years (Russell et al., 2001). The 
strong axisymmetry of the dipolar field of model J is in better agreement 
with the magnetic field of Saturn with a dipole tilt less than 1° (Russell and 
Dougherty, 2010). 

5.5 Convective motions within the conducting region 

In this work we have not taken into account the convective motions within 
the deep conducting region. Wicht et al. (2002) studied the linear stabil- 
ity of an imposed zonal flow in a spherical shell modeling the molecular 
hydrogen layer of Jupiter. They found that the critical Rossby number of 
the shear instability onset is almost independent of the Rayleigh number, 
which measures the strength of the convection. They concluded that the 
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shear instability is only weakly modified by the presence of convection. On 
the other hand, they showed that the convection onset is strongly influ- 
enced by the presence of the zonal circulation, with the convection either 
enhanced or damped depending on the direction of the shear. However, 
their study is linear, and so the results cannot be extrapolated beyond the 
weakly non- linear regime of convection. Whether or not our results apply in 
the presence of convection is a subject for future studies. In the presence of 
convection (which produces strong zonal motions and Rossby waves), and 
even for a convectively-driven dynamo (see for instance Aubert, 2005; Grote 
and Busse, 2001), the mechanism described here may still impose a simi- 
lar relationship between the magnetic field morphology and the zonal wind 
profile. 
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A Linear code used to compute the hydrodynam- 
ical instability threshold 

In order to study the linear stability threshold at very low Ekman numbers, 
we designed a linear code derived from Gillet et al. (2011). This three- 
dimensional spherical code uses second order finite differences in radius and 
pseudo-spectral spherical harmonic expansion. The linear perturbation u of 
the imposed background flow U is time-stepped from a random initial field 
with the following equation: 





— e z + V x U xu + UxVxu - Vp, 



(21) 



39 



together with the continuity equation V • u = 0, which allows us to eliminate 
the pressure term by using a poloidal-toroidal decomposition. The left hand 
side of equation 21 is treated with a semi-implicit Crank-Nicolson scheme, 
whereas the right hand side is treated as an explicit Adams-Bashforth term. 
Thanks to the cylindrical symmetry of the base flow U, all azimuthal modes 
m of the perturbation u are independent, and we can compute them sep- 
arately. The coupling with the background flow and the Coriolis force are 
handled in physical space, but a very fast implementation of the spherical 
harmonic transform (SHTns library) makes the code quite efficient. 

In order to determine the stability threshold at E = 10~ 7 , we used 
350 points in the radial direction, and spherical harmonics truncated at 
l-max = 300. We use no-slip boundary conditions. 
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